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ABSTRACT 

We test a particular theory of dark matter in which dark matter axions form ring 
“caustics” in the plane of the Milky Way against actual observations of Milky Way 
stars. According to this theory, cold, collisionless dark matter particles with angular 
momentum flow in and out of the Milky Way on sheets. These flows form caustic rings 
(at the positions of the rings, the density of the flow is formally infinite) at the locations 
of closest approach to the Galactic center. We show that the caustic ring dark matter 
theory reproduces a roughly logarithmic halo, with large perturbations near the rings. 

We show that the theory can reasonably match the observed rotation curve of the Milky 
Way. We explore the effects of the caustic rings on dwarf galaxy tidal disruption using 
N-body simulations. In particular, simulations of the Sagittarius dwarf galaxy tidal 
disruption in a caustic ring halo potential match observations of the trailing tidal tail 
as far as 90 kpc from the Galactic center; they do not, however, match the leading tidal 
tail. None of the caustic ring, NFW, or triaxial logarithmic halos fit all of the data. The 
source code for calculating the acceleration due to a caustic ring halo has been made 
publicly available in the NEMO Stellar Dynamics Toolbox and the Milkyway@home 
client repository. 

Subject headings: Galaxy: structure - Galaxy: kinematics and dynamics - Galaxy: 
stellar content 


Introduction 


1.1. The caustic ring dark matter halo 


T he caustic ring dark matter theory has been developed by Pierre Sikivie (ISikivie et al. 


1995 


Sikiviel 120031 : iNataraian fe Sikiviel 120071 : 1 Duffy &; Sikiviel 120081 : ISikiviel 1201 ll : iBanik &; Sikiviel 120131 ) 
ov er the past 20 years. A l thoug h Sik ivie started from the spherically symmetric self-similar models 
of lFillmore fc Goldreichl ()1984l ) and iBertschingen (j 19851 ). it is important to note that the Sikivie 
model has evolved to include physics of dark matter particles beyond a simple gravitational inter¬ 
action. This physics results in dark matter flows that rotate with the stellar disk. 
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The caustic ring model (Duffy & Sikivie 2008) is consistent with a Universe that expands 
with the Hubble flow, and has power law density fluctuations which create the seeds of individual 
galaxies. The dark matter starts out on these density fluctuations with temperature zero, so that 
the velocity dispersion at any point in space is zero. As the dark matter falls into a Galaxy it 
gains velocity, but since it initially occupies only three dimensions of six-dimensional phase space 
(three spatial dimensions and three velocity dimensions), it falls in on discrete flows that occupy 
only three dimensions of the six-dimensional phase space. The infalling flows, like the particles of 
which they are composed, fall in toward the galaxy center until they reach their closest approach, 
then fly out until all of the kinetic energy has been converted to gravitational potential energy, and 
then turn around and fall back in. 

For special locations in the Galaxy, the dark matter will form caustics, where the density is 
infinite in an infinitely small region of space (for non-zero velocity dispersion, the caustic will not 
be as sharply defined). These locations are the turnaround points on the orbits of the infalling 
sheets of dark matter. The caustics are rings at t he i nner turnaround radius, and spheres at the 
outer turnaround radius. In a self-similar model (Sikivie 20llj), the position of the caustics move 
outward with time, as initially more distant material collapses into the potential well of the galaxy. 
The number of caustics increases as more matter falls into the galaxy (with infall times of about 10 8 
years); each time the initial infalling particles reach perigalacticon, a new ring caustic forms near 
the center of the galaxy, at the original Galactocentric distance of that earliest accreting matter. 
Thus, the number of caustics grows with time (about one additional caustic for each 10 8 years); 
the most recently accreted matter is associated with the outermost ring. The locations of the ring 
caustics depend only on the slope of the evolved power spectrum of density perturbations on galaxy 
scales, the rotation speed of the galaxy, and the age of the Universe. For the Milky Way (rotation 
speed of 220 kms -1 at 8.5 kpc from the Galactic center) and the current standard cold dark matter 
(ACDM) model of the Universe (age = 13.7 Gyr), the inner ring caustics are predicted to be at 
radius: 

40 kpc 


n 


in = 1, 2,3,4,...) 


( 1 ) 


The first four caustics in the Milky Way are ring s at predicted distances of 40 kpc, 20 kpc, 13 kpc, 
and 10 kpc from the center of the Milky Way (Duffy &; Sikivie 20081). In cross section, these rings 
have a tricusp shape, as shown in Figure CD (top panel). 


1.2. Caustic rings and hierarchical structure formation 


Although the caustic ring model (we will also refer to it as the caustic ring halo) does not 
include the effects of hierarchical merging, it is plausibly consistent with the merger history of 
the Milky Way. Large galaxy merg ers w ill destr oy dark matter caustics similarly to the way 
galaxy mergers disrupt disks of stars (Toth & Qstrikerl 1992), and any other coherent structure of 
collisionless particles. However, it is very important to note that any dark matter that falls in 
after the collision will remain intact on the 3D sheets in 6D phase space. A merger produces a 
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local thickening in the ph ase spac e sheet , and a velocity dispersion associated with the flow of dark 
matter in a caustic ring. Sikivie (120031 ) estimates the spread of the caustic ring radius, 6a, due 
to the velocity dispersion, 6v, of the flow of dark matter. Minor mergers will cause the caustics 
to spread and become fuzzy, but the caustic will become sharp again after the merger event, as 
materia l from the cold flow fills the caustic in again. Minor mergers like the Sagittarius dwarf 
galaxy (|lbata et ah 1 99 4) that is currently being disrupted by the Milky Way’s tidal forces would 
not affect the structure of the caustics. 

The relationship between the redshift (z) at which the caustic ring forms and ring number (to) 
is approximately (P. Sikivie, private communication): 


z(n) = (3.8to - 0.9) u ' 4 - 1 


(2) 


Therefore the last (to = 1 is the last ring to form) fourteen caustic rings (to = 1,2,....14) were 
formed after z = 3.9, and the last two caustic rings (to = 1,2) were formed after z = 1.1. If there 
were no major mergers since 2 = 1.1, then we would expect to see a ring caustic (to = 2) at 20 kpc 
from the Galactic center, just where we see the Monoceros Ring of stars. 

In the Milky Way in particular, there is no evidence for a major merger in recent history. In 
particular, the boxy bulge of t he M ilky Way is cons istent with formation from disk instabilities 
rather than from mergers ( Ness et al.ll2013 : Shen et al. 2010 1. Even within the ACDM mo del, galax¬ 


ies w it h the lowest ratios of bulg e to total luminosit y have no major mergers at z < 2 (IMartig et al 


20121 s ) . iGilmore et al.l ([20021 ) and IWvse et al.1 (|2006l l have suggested that the thick disk may be the 
result of the last major merger of the Milky Way, 10 — 12 Gyr ago (1.8 < z < 4). If this is true, the 
four ring caustics, at 40 kpc, 20 kpc, 13 kpc, and 10 kpc from the center of the Milky Way, should 
be undisturbed. 

The details of the caustic ring model are still evolving; some of the ef fects of the B ose-E instein 
condensate (BEC) model are still being incorporated in the model (jBanik &, Sikivie 2013). The 
axions in this special degenerate condensed state (the de Broglie wavelength, %/p, must be large 
compared to inter-particle dimensions) set up a vorticity that is more permanent than the original 
model, when the axion self-interactions are included (the particles share angular momentum among 
each other). Just as water going down the drain of a sink has a curl (quasi-rigid rotation with the 
inner particles rotating faster than the outer particles), axions in a BEC also have a curl. In contrast, 
the velocity field for WIMPs has zero curl. Although the original density fluctuations are spherical 
by construction, the vorticity of the axions will result in an oblate dark matter halo. In addition, 
the baryons and the BEC dark matter may thermalize, thus explaining why disks have such a high 
angular momentum. Although very recent ACDM cosmological simulations can reproduce disks 
that are plausibly like those observed, it is unclear what p rocess gives t he disks rotation, while 
leaving the dark matter halo without significant rotation (see lMarinacci et al.ll20 14l. In this paper , 
we do not include the modifications to the caustic ring model suggested bv lBanik fe Sikiviel ([20 131) 
because their implications on the Galactic potential have not yet been worked out. 
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1.3. Implications of a caustic ring halo on the nature of dark matter 


The density fluctuations that are the precursors of galaxies in the standard caustic ring halo 
are by construction spherical. If the dark matter is additionally without rotation, then the dark 
matter falls in along radial orbits, and the inner ca ustic forms a cusp a t the center of each galaxy, 
with a density that diverges as 1/r 2 at the center ( Sikivie et al. 1 99a ). It is only when the dark 
matter has net rotation that the caustics form tricusp rings in the plane perpendicular to the axis 
of rotation, at the positions of the inner turnaround radii of the dark matter flows. This is in sharp 
contrast to ACDM simulations, in which the angular momentum of the dark matter distribution 
in most simulated galaxies is small. It requires a special form of dark matter to form galaxies with 
sufficient angular momentum to form the proposed caustic ring halos. 


Sikivie (J201 lj) showed that the dark matter distribution in galaxies has sufficient angular 


momentum if the dark matter is composed of axions rather than WIMPs or sterile neutrinos. The 
difference arises from the fact that the axions form a Bose-Einstein condensate (BEC) in which 
the particles exchange angular momentum through gravity. In general, the angular momentum of 
a galaxy as a whole is thought to be caused by a tidal torque from nearby protogala xies ea rly in 
its formation , before the galaxy has condensed out o f the expansion of the Universe (IHovld 1 195.31 : 
Peebles!Il969l : iPeebles!ll9Tll : lEfstathiou &; Joneslll980l gives a review!, though the dark mat ter can 
pick up net rotation through hierarchical merging (jYitvitska et al.ll2002l ; iDeason et al.ll201ll ). Since 
dark matter is typically modeled as collisionless particles, one would not normally expect that the 
dark matter would be susceptible to a tidal torque. However, axions in a rethermalizing Bose- 
Einstein condensate will react to a tidal torque by going to the lowest energy state compatible with 
the acquired angular momentum. That state is a state of rigid rotation on the turnaround sphere. 
The su r prising su g gest ion that axions would form a Bose-Einstein condensate was confirmed by 
Saikawa & Yqmaguchi (2013). If solid evidence of large angular momentum in the dark matter 
component of the galaxy is obtained, that would be an argument in favor of axion dark matter, or 
at least for additional dark matter physics. 


1.4. Evidence for a caustic ring dark matter halo 


Some observational evidence that caustic rings might exist in galaxies has been put forward by 
previous author s. This evidence include s: the statistical di stribution of bumps in the rotatio n curves 
of 32 galaxies (IKinnev fc Sikivie! [2000, using data from iBegeman et al.l Il99ll. and fr om I Sanders 
199 (tI); the distribution of bum ps in the rotation curve of the Milky Way (ISikiviel 120031 . using data 
from Oiling Rr. Merrifieldi 2000); a triangular feature in the Infrared Astronomical Satellite (IRAS; 
http://skyview.gsfc.nasa.gov) map in a direction tangen t to the nearest (n = 5) caustic ring, and 
aligned with the expected position of th e tricusp caustic (Sikivi3 2003 ): a nd the Monoceros Rin g at 
the position of the n = 2 caustic ring ( Nataraian <fe Sikiviel I2OO7I . citing Newberg et al. 2002 ). In 
order to match the rotation curves of the 32 external galaxies with the caustic ring model (which 
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has only been worked out for the Milky Way rotation speed), the rotation curves were rescaled by 
(220 km s _1 ) v, where v is the actual rotation speed of the galaxy. After this rescaling, there 
were 3 sigma detections of discrepant velocities in the rotation curve at (rescaled) radii of about 20 
and 40 kpc from the galaxy centers. The need for rings of da rk matt e r to e xplai n the in ner part of 
the M ilky Way galaxy rotation curve was also noted by ISofue et al.1 (|2009f ) and Ide Boer fe Weber 
(2011), who proposed two rings of dark matter at distances of 4.2 and 12.4 kpc from the Galactic 
center to explain dips in the rotation curve. The rings are close to Sikivie’s prediction of an n = 3 
caustic ring at 13.6 kpc and an n = 10 caustic ring at 4.3 kpc from the Galactic center, although 
it seems reasonable that the n > 10 rings might not be visible in the Milky Way rotation curve 
(because the distance between rings is less than a kpc or they were destroyed by ancient mergers). 
Possibly, the caustic ring model could be improved; for instance it does not currently include the 
presenc e of a bar in the disk or the effects of hierarchical merging. It is also interesting to note that 


Sikiviel (1200311. using data from the Massachusetts-Stony Brook North Galactic Plane CO survey 


( Clemens 


1985; 


Sanders et all Il986l : IClemens et al.1 1 19861 ). found evidence for about ten ripples 


in the rotation curve between 3 and 8.5 kpc from the Galactic center, which is about what was 
expected from caustics. The triangular feature in the IRAS maps, if caused by a tricusp caustic, 
implies p ~ 130 pc and q ~ 200 pc, where p and q are the horizontal and vertical sizes of the 
caustic. 


2. Calculating the Caustic Ring Gravitational Field 


Note that th ere is no way to add the caustic rings to a standard halo density model such as an 
NFW profile (INavarr o et al. 1996|), because the acceleration calculated in the caustic ring model is 
not the acceleration due to an excess of matter on a ring; it is rather the acceleration due to one 
sheet of dark matter falling in to the Milky Way. To derive the gravitational held for the entire 
halo, one must perform a vector sum of all of the accelerations due to all of the 20 (n = 1, 2,3,4,...) 
sheets of matter that have fallen in within the last 2 x 10 9 Gyr (previous infalls are expected to be 
a small fraction of the total halo mass). All of the halo dark matter, including the ring caustics, 
is included in this final calculation of the acceleration. There is no freedom to adjust the overall 
shape of the halo; it simply falls out of the sum of the integrals. 

Since the caustic ring model is axially symmetric, the equations for dark matter density and 
gravitational held in the Galaxy are conveniently expressed in cylindrical coordinates (p, z ), centered 
on the Galactic center. There is one equation for the density at a position near a caustic ring and 
another for the density at a position fa r away from a causti c ring. The equation for the density 
near a caustic is given by Equation 2 in lNatarajan Sz Sikivie (2007): 

N(p,z) 


d(p,z) 


1 


£ 


p 


dM 

dQdt 


;a,r) 


cos a 


\D(a,T)\ 


(a j(p,z),Tj{p,z)) 


( 3 ) 


Here the sum is over the number of hows of axions at (p, z), jctm ( a > r ) is the rate at which mass ( dM ) 
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falls in to the position (p, z) per unit solid angle and unit time (dQdt), and D(a, r) is the Jacobian 
determinant of the transformation between the coordinates (a,r) and ( p,z ). Each of t hese quan ¬ 
ti ties is a complex function of constants and coordinates in (a,r), defined in Natar aian fc Sikivie 
(120071 ). a is an angle associated with the particle’s last turnaround and r is the time when the 
particle crosses the z = 0 plane. The density of dark matter near a caustic ring is very large 
because there are many particles at turnaround points on their orbits around the Galaxy. The near 
equation shows this behavior since the density d(p, z ) approaches infinity when D(a, r) approaches 
0. When D(a,r) = 0, the transformation between the coordinates (a, r) and (p,z) gives a set of 
points (p,z) outlined by a tricusp shape. The tricusp shape is the blue solid line in Figure [[] (top 
panel). 


The equation for the density far from a caustic is given by Equation 3.45 in Duffy & Sikivie 


(2008)): 


1 dM 

d(p,z) far = -- 


(4) 

v dVldt y 7 (r 2 — a 2 ) 2 + 4 a 2 z 2 ’ 

where r = \J p 2 + z 2 , v is the speed of the flow near the caustic, ^is the mass infall rate per unit 
solid angle and unit time, and a is the caustic ring radius from the Galactic center. The equations 
given above for density apply to one caustic ring and its associated infalling and outgoing axions. 
The predicted radii (a n ), flow speeds (v n ), and infall rates ( 4^1 p .), for the first 20 caustics in the 


disk of the Galaxy are listed in Table 3 of Duffy & Sikivie |2008). These parameters will now be 
used with the subscript n that specifies the flow number (n = 1 to 20). 

The gravitational field for each caustic ring is similarly expressed by a near and far equation. 
If the curved nature is neglected near the ring, it can be modeled as a straight tube. This leads to 
a gravitational field of the following form: 


g(p, z) = —2 G J dp'dz' d(p', z') 


>- p')f>+ (z ~ z')z 


(p - p') 2 + {z- z') 2 ’ ^ 

( Nataraian fc Sikiviel 2007 1 where G is the gravitational constant and d(p',z ') is the dark matter 
density (Equation [3]). This equation describes the gravitational field fro m the fl ow o f all axions 


flowing into and out of the ring, and not just the ring structure itself. Tam (2012) presents a 


method to calculate this rather difficult integral anal ytical ly. After transforming coordinates and 
substituting d(p,z ) from Equation [3] into Equation [Sj Tam ( 201 21 reduces g into the form: 


5(P-z) 


8t tG dM 
pv n dVLdt 


[I P (p,z)p + h(p,z)z\ 


( 6 ) 


where IJp, z) and I z (p,z ) are integrals defined in Tam ( 20121 ). Here we replace the b in Equation 
16 from T arnl (120121 with v„. b is of the same order of magnitude as the speed of the particles in 
the caustic ([Nata raian fc Sik ivie 20071 1. Figure Q] (top panel) shows the vector representation of 
g(p, 2 )near for the n=2 caustic ring which lies at the radius a n of 20.1 kpc. The field points towards 
the caustic in all directions. The tricusp is symmetric around a central point where the field goes 
to zero. For the n=2 ring (Figured) top panel), the central point lies at (p,z) = (20.35,0) kpc. 
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Fig. 1.— (Top) Vector representation of gravitational field g nea r for the flow of dark matter near 
the n = 2 caustic ring, p and z are Galactocentric cylindrical coordinates. The blue solid line 
outlines the ‘tricusp’ shaped ring. (Bottom) Vector representation of (?haio for the sum of the flows 
from all caustic rings (n = 1 — 20) in the Milky Way. The vertical blue lines are at the postions of 
the caustic rings, p and 2 are Galactocentric cylindrical coordinates. 
























The gravitational field at a position far from the caustic is: 


g(p, *)far = fo 9 ", x [( r 2 - a 2 n )pp + (r 2 + a 2 n )zz \, 
s(2a 2 n + s) 

_ 4vr G dM 
' " ~ v n dfldt 

where r = \J p 2 + z 2 and s = (r 2 — ) 2 + 4a^z 2 (P. Sikivie, private communication). 


( 7 ) 

( 8 ) 


The caustic ring halo acceleration vs. Galactocentric distance in Figures 0 (bottom panel) and 
[2] was calculated as follows: if the distance (p, z) is inside the tricusp ring, use the field near the 
ring (Equation [6]); if the distance (p, z) is outside the ring use the field far from the ring (Equation 
0 . The gravitational field at a distance (p, z) is the vector sum from all caustic rings (n = 1 to 20): 


20 

g(P, -2)halo = 5n( P , Z). (9) 

72 — 1 


We spent considerable effort finding the spatial limits that made the near-caustic function match 
seamlessly to the far -caustic function, and finding and deciphering the published model parameters 
for the Milky Way (jDuffv Sz Sikivid 2008i ). There were essentially no tunable parameters in the 
model. Figure 0 (bottom panel) shows the gravitational field vectors for ghalo as a function of 
position in the Milky Way. The blue lines show the positions of one cross section of the caustic 
rings in the pz plane. The field generally points towards the Galactic center except near the 
innermost ring, where the field points away from the center. This is because there is no mass inside 
the n = 20 ring in our model, so objects not exactly at the galaxy center are pulled towards the 
innermost rings. The contribution to the halo mass from the additional caustics for n > 20 is 
assumed to be negligible. These additional caustics have likely been rearranged by early merger 
events. Although these figures show the cross section of the rings in the pz plane, the tricusp is 
barely resolved on this scale. This is because the average horizontal size of the tricusp is on the 
order of ~ 0.4 kpc. Note the vertical size of the tricusp has been exaggerated in Figure 0 (blue lines 
in bottom pane l); this average s i ze is ^ 0.5 kpc. The sizes of the tricusp are not predicted by the 
caustic theory. Duffy_& Sikiv ie (20081 ) estimate these sizes from rises in the Milky Way rotation 
curves, which will be discussed later in this paper. The reasoning is that (at z = 0) there are 
discontinuities in the rotation curve at positions that lie on the boundary of the ring. For example, 
in Figure 0 the tricusp boundary for z = 0 lies at p = 20.1 kpc and p = 21.1 kpc. This implies 
there will be discontinuities in the rotation curve at p = 20.1 kpc and p = 21.1 kpc. 


The source code for calculating the caustic ring halo acceleration given a position has been made 
available on the Milkyway@home GitHub repository (http: //github. com/Milkyway-at-home/milkywayathome_ 
The caustic ring halo has also been released as a gravitational potential module in the NEMO stel¬ 
lar dynamics toolbox (http://bima.astro.umd.edu/nemo/). Orbits and N-body simulations can 
be run on particles in a caustic ring halo potential using both NEMO and Milky way @home. 
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3. Comparison of the Caustic Ring Halo with Other Halo Models 


The halo potential models (NFW and logarithmic) that will be compared to the caustic ring 
model throughout this paper are listed below. X, Y, and Z are Galactocentric Cartesian coordinates 
for all models used in this paper. 


T he equation for the spheric al NFW model (jNavarro et alJll99a ) was taken from Equation 
2.67 in Binne y X Tr emaine (2008): 


^nfw = -47rGp 0 r 


2 ln( l + r/r s ) 
r/r s 


( 10 ) 


where r = \JX 2 + Y 2 + Z 2 . po is a characteristic overdensity of the halo and r, is the scale radius. 
The parameters we used for the NFW model were adopted from Xue et ah (j2008l ). In t his paper 


t he o bserved radial velocities of BHB stars from the Sloan Digital Sky Survey (SDSS; York et al 


20001) were compare d to those from simulations of the for mat ion of galax ies like the Milky Way. 


In Xue et al. (120081 ). po is defined as p s (Equation 14 in Xue et al. 20081). and r s = r v i r /c. The 
characteristic density p s depends on the critical density of the Universe, the matter contribution 
to the critical density, and the critical overdensity at virialization. r v i r is the virial radius and c is 
a dimensionless concen trat ion param eter. A Milky Way model including an NFW halo was fit to 
the rotation curve. Xue et al. (2008) found fit values for the virial radius and concentration. We 
use those fit parameters here for the NFW model. The values for the parameters in Equation [TO] 
are listed in Table [T] 


The parameters for a logarithmic (log) halo were taken from Law et ah ( 200a ). Law et al 


(120 0a) used a log halo in the Galaxy model to fit the orbit of the Sagittarius dwarf galaxy tidal 
stream. We discuss the Sagittarius stream in Section 5. Uhaio is the halo speed, q is the Z-direction 
flattening, and d is the halo scale length. The functional form of the logarithmic halo is: 


3*iog — 0 2 halo In (x 2 + T 2 +-p- + d 2 ^j . 


( 11 ) 


The parameters for a triaxial halo were taken from Law Sz Maiewskil (2010). Uh a io> T an< l ^ are 
defined the same as they are for the logarithmic halo, and the functional form is: 


<F 


triaxial — v t, halo 


In 


7 2 

C x X 2 + C 2 Y 2 + C 3 XY + - 2 - 

9t 


+ d‘ 


( 12 ) 


The constants C\, C 2 , and C 3 are defined in Equations 4-6 of Law fc Majewski (j201y). These allow 
for rotation ab out t h e halo Z axis and flattening in the X and Y directions. These parameters 
were also fit bv ILaw Sz Maiewsk i Q2010|) to match the orbit of the Sagittarius tidal stream and are 
listed in Tabled] 


The gravitational acceleration of the NFW and spherical (q\ = 1.0) log halo models in Figure 
[2] was calculated using g = — V4>. It is important to note that unlike the NFW, log, and triaxial 
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Fig. 2.— Gravitational acceleration vs. distance from the Galactic center. The cyan line shows 
the acceleration in the p direction for points in the Galactic plane {z = 0) as a function of p. Note 
the bumps at positions where caustic rings are predicted. The bumps near the first three caustic 
rings at 40.1 (n = 1), 20.1 (n = 2), and 13.6 (n = 3) kpc from the Galactic center are labeled. The 
red line shows the acceleration in the z direction for points along the axis of symmetry (p = 0) as 
a function of z. Since there are no caustics along the axis of symmetry of the Galaxy, no peaks 
appear in the red curve. The acceleration of the NFW (blue) and logarithmic halo (green) are 
shown for comparison; since these two potentials are spherically symmetric, the acceleration as a 
function of p and z are the same. 
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models, the caustic ring model parameters are predicted directly from theory. The caustic ring 
model parameters a n , v n , and \ n (Equations I61IH1) follow from solutions to the equations of 
motion for axion flows falling into the Galaxy. The initial conditions are specified by the initial 
dark matter mass distribution, the rotation speed of the Galaxy, and the initial dark matter flow 
angular momentum. The ‘tunab le’ parameters a re the angular momentum parameter jmax and the 
size p of the caustic ring (Duffy & Sikivie 2 0081) . T he angular m omentum parameter j ma x f° r the 
Milky Way was fit t o the rotation curve b y lSikivie et al.1 ill9951 ). The model predicts the radii for 
the first 20 caustics (Duffy & Sikivie 20081). The contribution to the halo mass from the additional 
caustics for n > 20 is assumed to be negligible. The si zes p o f t he caus tic rings are inferred from 
the analysis of bumps in rotation curves. Table 4 of Duffy fc Sikiv ie (120081 ) lists the positions 
of observed rises in Milky Way rotation curves corresponding to the positions of the n = 3 and 
n = 5 — 14 caustic rings. There should be an upward kink in the rotation at the position of one 
edge of the caustic and a downward kink at the position of the other edge of the caustic. In Figure 
Q] (top panel), this corresponds to the left edge at p = 20.1 kpc and the right edge at p = 20.4 kpc. 
For our c alculations of the caust ic ring model gravitational acceleration we use the p values from 
Table 4 of Duffy fc Sikivie ( 2008 ) for the n = 3 and n = 5 — 14 rings. For n = 1 — 2, 4, and 15 — 20 
we use Pn/o-n = 0.05 which is the average value of the ratio p n /a n for the other rings. 


The caustic ring halo acceleration is compared to that of an NFW halo and a spherical log halo 
in Figure [2j The acceleration for the NFW and log halos are given as a function of distance from 
the Galactic center. For the caustic ring halo, which is not spherically symmetric but is axially 
symmetric, we show acceleration in the plane (z = 0) as a function of cylindrical radius p, and along 
the axis of symmetry (p = 0) as a function of z. In either case the direction of the acceleration 
is towards the Galactic center. The obvious feature is that in the Galactic plane there are large 
jumps in the acceleration at the positions of the caustic rings. Since the NFW and log halos are 
spherically symmetric, there is no difference in these curves for z = 0 and p = 0. 


Given the complexity of calculating the acceleration, it is extremely gratifying and reassuring to 
discover that the resulting potential was plausibly similar to familiar Galactic potentials. Note that 
the caustic ring halo, overall, is reasonably similar to a log halo, but with jumps in the acceleration 
near the caustics. The shape of the dark matter profile differs tremendously from the NFW profile 
near the Galactic center; there is no strong central cusp in the caustic ring halo. In fact, it is just 
the opposite; there is essentially no dark matter in the center of a caustic ring halo. The sharp 
spike in Figure [2] near the center indicates an acceleration away from the Galactic center, towards 
the caustic rings. Figure Q] (bottom panel) shows the gravitational field vectors as a function of 
position in the Milky Way, showing the positive acceleration along the radius inside the innermost 
ring caustic. This affect disappears when bulge and disk potentials are included, as we shall see 
later. In the real galaxy, we can imagine that the caustic ring halo is slightly blurred due to the 
non-zero velocity dispersion of the dark matter. Also, past major mergers would disrupt the caustic 
pattern near the Galactic center. 
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4. Comparison of the Caustic Ring Halo with the Milky Way Rotation Curve 


The potential of the Milky Way galaxy includes contributions from the dark matter halo, the 
disk of baryons, and the stellar bulge. In this paper, we use analyt ical di sk and b ulge potentials 
to model the baryons in stars and gas in the Milky Way. We use a lMivamoto &; N avail (|1975l l disk 
and Hern auisti (1990!) bulge: 


<f>disk = - 


GM disk 


Y^X 2 + Y 2 + (a + y/Z 2 + b^y -’ 


(13) 


^bulge — 


G'-^buIge 
r + C 


(14) 


where r = \/X 2 + Y 2 + Z 2 . i s k and Mb u i ge are the masses of the disk and bulge, a and b are the 
disk scale length and height, c is the bulge scale radius. X, Y, and Z are Galactocentric Cartesian 
coordinates for all models used in this paper. We will compare the results of three cases for the 
fixed analytical potential of the Milky Way: 


1- *1* Galaxy — ^jjernquist bulge T ^Miyamoto—Nagai disk Y Realistic halo 

2. d'Qgjaxy 'hjjernquist bulge Y Y Miyamcto—Nagai disk Y ^logarithmic halo 

3. Y(h|] ax y — 4?hlernquist bulge Y ^Miyamoto—Nagai disk Y Ytriaxial halo- 


The three halo models used here are discussed in Section 3. The rotation curve was obtained as 
follows. For circular orbits the centripetal acceleration at a distance r from the Galactic center 
is g{r) = V ' r J , the square of the velocity at at r divided by r. We get the acceleration from 
from the Miyamoto-Nagai disk, Hernquist bulge, log halo, and NFW halo potential models using 
g(r ) = — V<h. For the caustic ring halo, we already have the acceleration. To get the con trib ution 
from the gas in the disk, we used the velocity points of the gas disk model from Oilin g fc Merrifield 
(120001 ). Our g gas is a fifth order polynomial fit to the HI gas model and a sixth order polynomial 
fit to the Hu gas model. Since we have the acceleration for all of our models, we can solve for v(r). 
We add all the contributions to the Galaxy acceleration gtotai = <?disk Y <?buige Y g ga s Y 5haio and plot 
v = y/Stotal • r for the rotation curve. 

We use a reduced y 2 as a measure of fitness of the rotation curve data points to the models: 


xl = 


1 


N — n — 1 


E 


model,i Vdata,i 

&v,i 


(15) 


where v is the velocity, a lKl is the error in the measurement of the data point, N is the number 
of data points, and n is the number of parameters being fit. We performed the fit by finding the 
optimum y 2 value for rotation curves produced with a given set of parameters for the disk, bulge, 
and halo models. Reduced y 2 values around 1.0 represent a reasonable fit to the data, implying 
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r (kpc) 


Fig. 3.— Rotation curve of the Milky Way. The sum of a caustic ring halo (lower solid line), 
Hernquist bulge (dotted), Miyamoto-Nagai disk (dashed), HI gas layer (dashed-dotted line), and H 2 
gas layer (dashed-dotted line) produce a rot ation curve (upper solid line) that matches observations. 
The blue data points are from ISofuel (|2012f ) . The d ashed-dotted line is a poly nomial fit to the gas 
model (open circles for H 2 and crosses for HI) from [Oilin g & Merrifieldl (120011 ). 
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most of the model point values are within the errors of the data point values. The parameters for 
the log (para meter valu es ado pted from Law et ah 2005, see Section 3), NFW (parameter values 
adopted from IXue et al. 2008, see Section 3), and caustic ring halo were not fit. 


In this case, the x 2 value now depends only on the disk and bulge parameters: x 2 (Md; s k, a, b, 
Mbuigej c). The fit parameters were found by parameter sweeps; varying one parameter while 
leaving the others constant and finding the chi-squared values. This consequentially finds the x 2 
as a function of one parameter. For example, if we have an arbitrary parameter p , we can find 
X 2 (p). If the plot of x 2 (p) is close to a horizontal line, the fit doesn’t depend on this parameter. 
Parameters that affect the fit produce x 2 {p) plots that are close to parabolas with a dip at the 
minimum chi-squared value. We found that the fit here depends on the disk mass (M^isk), disk 
scale length (a), bulge mass (Mbuige), and bulge scale length (c) but doesn’t depend on the value 
for the disk scale height ( b ). We have x 2 (Afdisk> a, Mb u i ge , c) with 4 free parameters. The number 
of data points is N = 121 and the number of parameters is n = 4 (Equation 1151) . The fit values we 
found for the disk and bulge parameters are shown in Table [2j 


In Figure 01 the sum of a caustic ring halo (lower solid line), Hernquist bulge (dotted), 
Miyamoto-Nagai disk (dashed), HI gas layer (dashed-dotted line), and H 2 gas layer (dashed-dotted 
line) prod uce a rotati on curve (upper solid line that matches observations. The blue data points 
are from Sofue (.2012). The data from Some (120121 ) is a compilation of previous rotation curve 


data from the literature with the velocities scaled to assume a rotation speed of 200 km s _i at 
8.0 kpc. The dashed- dotted li ne is a p olynomial fit to the gas model (open circles for H 2 and 
crosses for HI) from Oiling fe Merrifieldl (2000) as discussed above. The caustic ring halo rotation 
curve maintains kinks in the curve at the position of the caustics, similar to the acceleration plot 
(Figure [2]). Although the data and the model are similar, the observations are not detailed enough 
particularly for r > 10 kpc to correlate bumps in the rotation cu rve with rises in the data points, 
as was suggested as evidence for the existence of caustic rings by [Sik ivie (20031 ). 


Figure 0] shows rotation curves generated with the same disk, bulge, and gas components as 
Figure 01 except with different halo models. The rotation curve fit for 0 < r < 2 depends primarily 
on the contribution of the bulge model. The bulge mass and scale radius determines the fit to the 
region near the central peak at r ~ 0.3. The fit for 2 < r < 15 depends mostly on the contributions 
of the disk and halo. The caustic ring model has a negligible amount of mass within 2 kpc of the 
Galactic center and also does not contribute very much to the rotation curve for the inner Galaxy. 
This is shown in Figure 0] (top panel) where the disk curve is much higher than the halo curve for 
0 < r < 15. The same is true for the case of a log halo model (middle panel). The NFW halo 
conversely predicts a much larger amount of mass for the inner Galaxy. In Figure 0J the NFW halo 
curve is almost level with the disk model curve. This exposes a degeneracy in the choice of the disk 
or halo model to fit the rotation curve data. Here the choic e for a caustic r ing h alo or log halo are 
consistent with the maximum disk hypothesis (discussed in IBinnev &; Tremaine! 2008l). This is the 
idea that the disk should contribute most (50-90% of the total radial force) to the rotation curve 
of the inner Galaxy, and assumes that the mass-to-light ratio of the disk doesn’t depend on the 
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Table 1. Halo Model Parameters 


Model 

Parameter 

Value 

NFW 

Characteristic density 

p 0 = 4.197 x 10” 3 Mq pc -3 


Scale radius 

r s = 22.25 kpc 

logarithmic 

Halo speed 

id, halo = 114.0 km s -1 


Z direction flattening 

qi = 1.0 

triaxial 

Halo speed 

id, halo = 121.9 km s' 1 


Z direction flattening 

qt = 1-36 


Scale length 

d = 12.0 kpc 


Constant 

Ci = 0.992947 


Constant 

C 2 = 0.532153 


Constant 

C 3 = 0.114889 


Table 2. Rotation Curve Fit Parameters 




d*bulge T ^diskV 

^bulge T ^diskT 

*^bulge T ^diskV 

Parameter 


^caustic ring halo 

<hlog halo 

^NFW halo 

disk mass 

Miisk (M e ) 

9.89 x 10 10 

9.11 x 10 10 

7.33 X 10 10 

disk scale length 

a (kpc) 

5.0 

5.0 

5.0 

disk scale height 

b (kpc) 

0.26 

0.26 

0.26 

bulge mass 

Vfbulge (Mq) 

1.89 x 10 10 

1.89 x 10 10 

1.89 x 10 10 

bulge scale radius 

c (kpc) 

0.33 

0.33 

0.33 


xl 

1.04 

1.01 

1.19 
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radius. 


The fit values for the disk and bulge parameters are shown in Table [2j The main difference 
in the fits for the three panels in Figure [4] is the value of the disk mass (Mdi s k): 9.89xl0 10 , 
9.11xl0 10 , and 7.33xl0 10 M 0 . In Table El a, b, Mb u i ge , and c are identical for each column. 
Estimates for the Milky Way disk and bulge m asses vary dependi ng on which particular model and 
technique is used to fit the observational data. Dehnen & Binncy (Jl998|) used a Galaxy model with 
a spheroidal bulge, exponential disk, and spheroidal halo to f it values of Mdi sk ~ 5.13 x 10 10 M 0 


and Mbuige ~ 0.518 x 10 10 M 0 for the disk and bulge. Wi dro w _& Dubinski ( 200 51) fit values 
of M disk ~ 4.58 x 10 10 M 0 and Mb u i ge ~ 1.1 x 10 10 M 0 for a Galaxy model with a central 


black hole, Hernquist bulge, exponential disk, and NFW halo. McMillanl (2011) fit a value of 
Mdisk ~ 5.76 x IQ 10 Mm assuming a bulge mass of Mb u i ge = 0.897 x 10 10 for a Galaxy mo del 
with a lBissantz fc Gerhardl (2002) bulge, exponential disk, and NFW halo. iKafle et al.l 1 20141 ) fit 
values of M d i sk ~ 9.5^3 q x 10 10 M 0 and Mb u i gc ~ 0.91^Q3g x 10 10 M 0 for a Galaxy model with 
a spheroidal bulge, Miyamoto-Nagai disk, stellar halo, and NFW dark halo. The disk masses we 
found here are larger than what was found for an exponential disk but within the errors of what 
was found for a Miyamoto-Nagai disk. Our bulge masses are sightly larger than what was found in 
the sources listed here. 


5. The Sagittarius Dwarf Tidal Stream 


In this context, tidal streams are collections of stars that have been gravitationally stripped 
from dwarf galaxies that in turn are orbiting much larger galaxies. As a dwarf galaxy is tidally 
disrupted, two streams or “tails” of stars are formed. The leading tail moves ahead of the orbit 
of the boun d dwar f galaxy and the trailing tail follows behind it. The Sagittarius (Sgr) dwarf 


galaxy ( 

Ibata et al. 

1994 

) tidal stream is the largest and best studied tidal stream in the Milky 

Way galaxy ( 

Maiewski et al. 2003, 2004: Vivas et al. 2005: Johnston et al. 2005: Law et al. 2005: 

Law & Maiewski 

201 o|b 


The Sgr stream has been used to measure the Milky Way potential and thus constrain the 
shape of the dark matter halo. The data for Sgr tidal debris has however been difficult to fit to 
halo models. There have been claims for spher ical, prolate, and oblat e ellipsoid models. Since the 
Sgr debris is approxi mately confined to a plane, iMaiewski et al.1 (120031 ) suggested that the potential 
should be spherical. Vivas et al-J (2005) compared radial velocities and distances from the leading 
tidal tails to numerical simulations for different halo models. They found that prolate and spherical 
models fit better but still fail to reproduce the entire angular span of the stream. The tilt of the 
tidal tails with respe ct to the Sg r orbit al plane, however, suggests that the potential is oblate 
Johnston et al.l 1 200 51 ). lLaw et al.l (120051) conve r sely f ound the prolate models best fit the leading 
tail velocities. More recently, [Law Sz Maiewskil (20KJ) found a reasonable fit to both tidal tails if 
th ey assume a triax ial halo, with both the major and minor axes in the plane of the Milky Way. 


Vera-Giro & Helmil (2013) expanded upon this idea to fit a halo potential that varies as a function 
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Fig. 4.— Rotation curve of the Milky Way. Reproduction of Figure[3]with the logarithmic (middle) 
and NFW (bottom) halo shown for comparison. All three models are consistent with the measured 
rotation curve, but with different disk masses. 
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of Galactocentric radius. These authors obtain an oblate fit within 10 kpc of the Galactic center 
and a fit to the triaxial model for outer radii. Due to the controversy in modeling the tidal tails 
of the Sgr stream, it is worth asking whether substructures such as caustic rings could affect the 
shape and velocity structure of this tidal debris stream. 

As a test of the viability of the caustic ring model, we determine whether Sgr stream observa¬ 
tions are consistent with orbits and N-body simulations in a caustic ring halo potential. In Section 
5.1, we discuss the orbit fitting technique used, and the results of comparing the model to the 
observational data. After the orbits are fit, N-body simulations are run with those orbits. We will 
compare the N-body models to the data in Section 5.2, using three cases for the fixed analytical 
potential of the Milky Way: 


1- ^Galaxy — ^pjernquist bulge T ^Miyamoto—Nagai disk T ‘hcaustic halo 

2- d’Galaxy = ^Hernquist bulge T ^Miyamoto—Nagai disk T ^NFW halo 

3. ^Galaxy ^Hernquist bulge T T Miyamoto_Nagai di sk T ^triaxial halo 

The galaxy models used here are those that have been fit to the Milky Way rotation curve in Section 
3. The disk, bulge, and halo parameters are from Table [T| and Table [2j For the case of the triaxial 
halo potential, we use the disk and bulge parameters from the log halo fit in Table Q] and Table [2j 


5.1. Orbit Fitting 


An orbit (line through phase space) of a test particle through a particular external potential is 
specified by initial conditions for position and velocity. We generat e orbits using the mkorbit and 
orbint tools from the NEMO Stellar Dynamics Toolbox (ITeubenl 19951) - The mkorbit tool generates 
an orbit hie given a specified analytical external potential and starting values for position and 
velocity. The orbint tool integrates the orbit through a specified number of timesteps. We use the 
orblist tool to list the position and velocity values calculated at each timestep. 

We can constrain the potential of the Milky Way by matching model orbits with observational 
data from stars in tidal streams. We used a gradie nt descent method to fit Sgr stream data to 
orbits using methods described in section 6 of Willett et ah (120091). The NEMO code generally uses 
Galactocentric Cartesian coordinates for input of positions and velocities (x, y, z, v x , v„ , y?). For 


convenience, we will use the Sgr heliocentric coordinate system defined in lMaiewski et al.1 (120031 ) to 
designate positions (d sun , Aq, Bq). In this section and the next we adopt a distance of R & = 8.0 
kpc from the Sun to Galactic center. Now the initial conditions for an orbit are specified by: 

(c^sun; Aq, Bq, V x , Vy , V z ). 











-19- 



Fig. 5.— Galactic standard of rest velocities vs. orbital longitude (Sgr heliocentric spherical coor¬ 
dinate system) for the Sgr dwarf galaxy tidal disruption. The fo rward (dotte d line) and backwards 
(solid line) best-fit orbit to the 2MASS M giant data from Maiewski et al. (2003) (black dots) is 
shown overplotted with their corresponding N-body simulations (red dots). The orbits and N-body 
fits for the caustic (top), NFW (middle), and triaxial (bottom) halo are shown for comparison. 
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A 0 (degrees) 


Fig. 6.— Orbital latitude vs. orbital longitude (Sgr heliocentric spherical coordinate system) for 
the Sgr dwarf galaxy tidal disruption. The forw ar d (dott ed line) and backwards (solid line) best-fit 
orbit to the 2MASS M giant data from iMaiewski et al.l (120031 ) (black dots) is shown overplotted 
with their corresponding N-body simulations (red dots). The orbits and N-body fits for the caustic 
(top), NFW (middle), and triaxial (bottom) halo are shown for comparison. 
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We use the reduced x 2 as measure of fitness of a model orbit to stream data points: 

/ \ 2 

2 = V ( ^gsr model ,i ^gsr data ,i 


'X.'Ucrar 




G i) 


Xb = 


E 


B model,* ^data,i 
&BA 


and 


x 2 = 


N — n — 1 


(xL r + Xb)- 


(16) 


(17) 


(18) 


where u gsr is the Galactic standard of rest velocity, a is the error in the measurement of the data 
point, N is the number of data points, and n is the number of parameters being fit. Since the orbit 
is a line through phase space, the best fit is one where the orbit line goes mostly closely through the 
set of stream data points. The fit is therefore determined by the starting positions and velocities 
and the choice of the external Galactic potential model. Here we can therefore compare orbit fits 
in different external potentials. 

Since the current coordinates of the Sgr dwarf galaxy are well known, we fix the starting 
points for the distance and orbital longitude and fit the starting points for the orbital latitude 
and velocity (d sun . A 0 , B & , v x , v y , v~) = (28.0 kpc, 0.126°, Bq, v x , v y , v~). We are fitting 4 
parameters. The orbit was fit to 2MASS M giant angular position/velocity data from Table 3 of 
Maiewski et al. ( 2004 ). Only M giants with d sun > 13 kpc and —7.8 < Bq < 7.8 were selected for 
our fits. These points were split into 10° bins in latidude A 0 , and points with velocities greater 
than 3 cj of the bin average were removed fr om the data se t. The error of each velocity point was 


taken to be <7„ gsr ,i = 6 kms 1 as stated in lMaiewski et al.1 (12 0041 ). The error in latitude, OB,i = 3 


is the standard deviation of Bq in our data set. The number of data points is IV = 80 and the 
number of parameters is n = 4. We completed orbit fits for the three cases listed above for the 
external Galactic potential, which highlights the effects of varying the choice of the halo model. 

The gradient descent method takes an initial set of parameters (d sun , A 0 , B & , v x \, v y \, v z \), 
generates an orbit, compares the velocities from this orbit to Sgr M giant data, and calculates a x 2 
value, Xoid- Next the velocities are adjusted to new values (d su n , A 0 , Bq, v x i , v y 2 , v Z 2 ) by moving 
in the direction of the x 2 gradient, a new orbit is generated and x 2 e w calculated. If x 2 ew i s l ess than 
X 2 M , the new set of parameters are kept. If x 2 ew is greater than x 2 id> the old set of parameters are 
kept. The process is repeated until Xnew = Xoid> signaling that a minimum has been reached. The 
best fit parameters, (d sun , A 0 , Bq, v x , v y , v z ) = (28.0 kpc, 0.126°, 5®, bestfit, v x , bestfit, %bestfit, 
v z bestfit) correspond with a m inimu m x 2 value . The errors in the fit parameters are estimated 
using the method described in I Cole et al. ([2008). We assume the shape of the x 2 surface near the 
minimum is Gaussian. The curvature of a surface is found using the Hessian or curvature matrix 
H, a matrix containing all possible second derivatives of a fu nction . W e calc ulate the derivatives 
for the Hessian matrix numerically using Equation 23 from Cole et al. (2008). The inverse of the 
Hessian matrix gives the error with respect to each best fit parameter: 

v^-Lh-W 


(19) 
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The diagonal elements of V are the variances and the off diagonal elements are the covariances. 
The variance is normalized by the number of data points N. cr* = s/Va/N gives the error estimate 
with respect to each best fit parameter. 


The starting velocities for the gradient descent were random values chosen from a range: 
(v x , v y , v z ) = (150 to 250, -40 to 20, 150 to 250). These values were chosen to ensure the gradient 
descent algorithm begins close the orbit of Sgr that is known from previous publications. The gradi¬ 
ent descent however will not necessarily end with values in these ranges. A disadvantage to using the 
gradient descent method is that if the starting point is very far from the global minimum, the algo¬ 
rithm may become trapped in a local minimum. We start with these values to ensure the global min¬ 
imum is found. Our best fits are chosen as the run with the lowest x 2 value after 10 separate runs. 
The results for the fits in the 3 cases listed above are: (B Qi bestfit, v x , bestfit, % bestfit, v z , bestfit) = 

1. (0.028° ± 0.7, 210.2 ± 1.1 km s" 1 , -30.8 ± 1.1 km s" 1 , 212.7 ± 0.7 km s" 1 ), x 2 = 52.2 
(caustic) 


2. (0.622° ± 0.7, 222.9 ± 1.4 km s -1 , -32.8 ± 1.4 km s" 1 , 208.7 ± 0.8 km s" 1 ), x 2 = 45.8 
(NFW) 

3. (0.441° ± 0.7, 224.6 ± 1.4 km s -1 , -29.4 ± 1.4 km s -1 , 198.7 ± 0.9 km s" 1 ), x 2 = 53.7 
(triaxial) 


Figures [5] and [6] show the best fit orbits to the M giant data (black points) for cases 1-3. The 
trailing data lies on the backwards orbit (solid line) and the leading tail data lies on the forwards 
orbit (dotted line). We have plotted the Galactic standard of rest velocity (u gsr = x-v x +y-v y +z-v z ) 
vs. orbital longitude Aq in Figure [H The angle along the tidal stream, Aq, is zer o at the positio n 


of the Sgr dwarf galaxy and defined to increase in the direction of trailing debris ( Maiewski et al 


20031 ). The x 2 values for each fit are on the order 50 because of the outliers in the data from 
Aq = 0° — 50° and the spread in range of the leading tail velocity points. The two data points on 
the edge of the leading tail are about 50-100 kins' 1 away from the orbit fit lines, yet from visual 
inspection of the plot the fits in the three panels are reasonably close to the data points. The 
trailing tail orbits in all three potentials are similar. In the leading tail, the orbits in NFW and 
triaxial halos favor upper data points. Yet the caustic ring halo orbit moves through the center of 
the points. Since the x 2 values for all three fits are similar, perhaps more data is needed in the 
trailing arm to determine which fit is favored. We have also plotted the orbital latitude, B@, vs. 
orbital longitude, Aq, in Figure [6l B@ ~ 0 is in the plane of the Sgr dwarf galaxy orbit. The orbit 
fits in each panel are all well confined to the Sgr plane. The N-body simulations (red points) in 
Figures [5] and [6] are discussed in the next section. 


We estimate the radial period of the best fit orbits by calculating the time difference in going 
from the first backwards orbit apocenter to the first forwards orbit apocenter. We found the radial 
periods to be 1.09 Gyr (caustic), 1.18 Gyr (NFW), and 0.92 Gyr (triaxial). These values are similar 
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to those found in Law et ah (2005)): 0.85 Gyr, 0. 
potentials, respectively. 


Gyr, and 0.87 Gyr for oblate, spherical, prolate 


5.2. N-body simulations of the Sagittarius dwarf tidal stream 


In this section, we will compare the results o f N-bodv si mu latio ns of the tidal disruption of 
the Sgr dwarf galax y. We us e the gyrfalcON tool (Dehnen 2000 . 2001 . 2002) in the NEMO Stellar 
Dynamics Toolbox ( Teuben 1995 1 to run the N-body simulations. Dwarf galaxies were modeled by 
a Plummer sphere (Plummer 191 1): 


d*P = — 


GMp 

\Jr 2 + r 2 0 


( 20 ) 


with initial mass Mp = Mg gr] o = 1.8 x 10 9 M 0 and scale length r D = 1.0 kpc with the number of 
bodies N = 10 4 . We assume the initial mass is the sum of the light and dark matter components of 
the dwarf galaxy. The caustic ring model does not specifically address the dark matter distribution 
in dwarf galaxies, but it is reasonable to assume that, like the Milky Way, they would not include a 


central cusp as in an NFW profile. We u 
value of Mso-r n ~ 6.4 x 10 8 Afn used by 

se Ms„ r n ~ 1.8 x 10 9 Mm, which is slightly larger than the 
Law & Maiewski (20101 but smaller than 10 10 — 10 11 Mm 

found by 

Niederste-Ostholt et al. 

(2010) and used by 

Purcell et al. 

(2011 

). The exact initial mass 


of the Sgr dwarf galaxy is not important to the results of our simulation. 


The Plummer model is placed at a point along the backwards orbit (in the past) and evolved 
to the position of the dwarf galaxy at the present day. The model is evolved in an external Galactic 
potential corresponding to the three cases listed in Section 5.1. The fixed parameters for the 
potential are also the same as listed above. The simulations in Figures [5]|8] were run for a time 
equal to three orbital periods in the respective potentials: 3.28 Gyr (caustic), 3.54 Gyr (NFW), 
and 2.77 Gyr (triaxial). 

In Figures [5E1 the Plummer model (red points) is evolved along the best fit orbits determined 
in the previous section. The panels show the results in the three halo potentials: caustic (top), 
NFW (middle), and triaxial (bottom) as a function of the angle along the tidal stream, Aq. Figure 
[5] shows w gsr vs. Aq. The caustic ring halo plot has much lower velocities at the end of the leading 
(dashed line) tidal tail, and one can even imagine that the two M giant stars are each from a 
different apparent branch of the N-body simulation. The N-body models in all three halos match 
the trailing tail data well. The caustic and NFW debris matches the leading tail data well, whereas 
the triaxial debris matches five of the data points and misses the sixth one. More data points 
however would be needed to determine which fit is preferred in the leading arm. Figure [6] shows 
Bq vs. Aq. The N-body models in all 3 potentials are well confined to the Sgr orbital plane. 

In Figures [7] and [8] we compare the N-body simulations to a set of more recent Sgr data that 
extends farther along the leading and trailing tidal tails, in the range 150 < Aq < 250. Figure [7] 
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Fig. 7.— Galactic standard of rest velocities vs. orbital longitude for N-body simulations of the 
Sgr dwarf galaxy tidal disruption. Figure [5] is repro duced wit h recen t data for the trailing (blue 
squares) and leading (green squares) tidal tail from Belokurov et al. (20141). The error bars are 
shown but most of these are smaller than the size of the squares. 
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Fig. 8.— Ysgr,gc vs. Xsgr,gc (Sgr Galactocentric spherical coordinate system) for N-body sim¬ 
ulations of the Sgr dwarf galaxy tidal disruption. The blue dots show the galaxy evolved in a 
caustic (top), NFW (middle), and triaxial (bottom) halo pote ntial. T he forw ard (dotted line) and 
bac kwards (solid line) orb it is also plotted. Recent data from iNewbv et al.l (J2013j) (cyan squares) 
and Newberg et al. (2.003a) (red squares) is also shown for comparison. Note both the caustic ring 
halo and NFW halo tidal debris extends as far as the 90 kpc data, but the triaxial halo fits the 
leading tidal tail. 
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shows the simul ations overp lotted with recent velocity data for giant stars detected photometrically 
from SDSS by Belokurov et ah (2014|). Here the caustic and NFW debris match the trailing data 
(blue points) but the triaxial debris does not. Conversely, the triaxial debris matches the leading 
tail data (green points), but the NFW and caustic debris doesn’t for the points from Aq = 200° — 
250°. Figure [8] is plotted using the Sgr Galactocentric Cartesian coordinate system defined in 


Maiewski et al.l (120031 ). The center of the Sgr Galactocentric system is the point in the Sgr dwarf 
galaxy plane that is closest to the Galactic center. Because the orbit of the Sgr dwarf galaxy 
is nearly polar, (X sgr . ]gc , Y sgI . jgc , Z sgI . !gc ) are approximately in the directions (X, -Z, Y), where X 
points from the Sun to the Galactic center, Y is in the direction of the Sun’s motion, and Z in the 


distance observations from 

Newberg et al. 

(2003) and 

Newbv et al. 

(2013). I 

vfewbere: et al. 

(2003) 

measured the distances of BHB stars detected photometrically from SDSS. 

Newbv et al. 

(2013) 


measured Sgr stream center distances by fitting the spatial density of F turnoff stars in SDSS using 
the statistical photometric parallax technique. The caustic and NFW N-body models extend as far 
as 90 kpc from the Galactic center, much further than the triaxial N-body model. This is plausible 
because the caustic ring halo potential well is shallow (see Figure [2]), and the bodies can maintain 
enough kinetic energy to fly out to further distances. The caustic and NFW debris however doesn’t 
match the data for the leading arm positions. The triaxial debris is the best fit to the leading arm. 
We find here that none of the models are able to match both the leading and trailing tail data 
simultaneously. 


The orbital precession of the stream is found by calculating the change in the angle Aq between 
the leading and trailing apocenters. Since the precession depends on gravitational potential, it is 
worth observing how it varies in the three potentials we compare in this section. We fit 2nd order 
polynomials to the N-body model debris and then found the angle between the apocenters. The 



the rotation curve gets steeper. We can assume the caustic ring halo potential is approximately 
spherical since the curve along the radial and z directions have the same scale and general shape 
(see Figure El). By inspection of Figure (H it makes sense that the precession angle in the NFW 
halo is larger than that in the caustic ring halo since the NFW rotation curve is steeper. 


In summary, we show that N-body simulations of the Sgr stream in a caustic ring halo match 
most of the recent data for stream velocities and distances except on the edge of the leading 
tail at Aq = 200° — 250°. Neither of the caustic, NFW, or triaxial models are able to match 
velocity/distance data from the leading and trailing tails simultaneously. Lastly, the caustic halo 
N-body model can match Sgr distance observations in the trailing tidal tail, including debris out 
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to 90 kpc. 


6. The Effect of a Dark Matter Caustic Ring on a Satellite Passing Through It 

We ran N-body simulations of dwarf galaxies represented by Plummer sphere models on orbits 
that pass close to caustic rings, to see whether the high density areas would have a strong effect on 
tidal disruption. We found that although there is a large discontinuity in the acceleration near a 
caustic, the caustic ring gravitational field does not qualitatively affect the results of the simulation. 
There are, however, differences in the shapes of the orbits and resulting tidal tails in a caustic ring 
halo potential compared to the spherical log halo potential. In Figures 1911111 a Plummer model 
(Equation 1201) dwarf galaxy with Mp = 10 8 M@, r 0 = 0.5 kpc, and number of bodies N = 10 4 , is 
evolved in two different Milky Way gravitational potentials: 

!■ ^Galaxy ^Hernquist bulge T d* Miyamoto—Nagai disk “b ^caustic halo 
2. ^Galaxy — ^Hernquist bulge “b ^Miyamoto—Nagai disk “b ^log halo■ 

The parameters for these models are listed in Table Q] and Table [2j 

In Figures [9] and flOl the Plummer sphere is evolved along different orbital paths that move 
near the n = 2 caustic ring at 20.1 kpc from the Galactic center (left boundary of tricusp in 
Figure H]). The simulations start along the backwards orbit (solid line) and evolve to move along 
the caustic ring and then move along the forward orbit (dotted line). For Figure [91 the orbit is 
designed to move in the +Y direction through the side of the caustic ring in the plane of the Galaxy 
at ( X , Y, Z , v x , v y , v z ) = (20.1 kpc, 0 kpc, 0 kpc, 0 km s _1 , 200.0 km s _1 , 0 km s _1 ). X, Y, 
and Z are Galactocentric Cartesian coordinates with the position of the Sun taken as R@ = 8.0 
kpc or ( X , Y, Z) = (—8.0, 0, 0) kpc. The simulation time is 2 Gyr. The orbit in a caustic ring 
halo (top left panel) is approximately circular while the log halo orbit (bottom left panel) has some 
precession. The dwarf galaxies also travel to different positions after 2.0 Gyr. When we look at 
histograms of the number of bodies along Galactic longitude (l is convenient since we are looking 
in the Galactic plane) we find that the dwarf less disrupted in a caustic halo. The central peak 
in a caustic halo has about 50% of the total number of bodies and the log halo has about 30% of 
the total. We also ran simulations where the Plummer sphere ended in the same position for the 
two potentials; those are not shown here. In this case we found the same behavior where the dwarf 
galaxy disrupts more slowly in a caustic potential. Here the dwarf galaxy may be on a stable orbit 
around the ring and the bodies are stripped from the galaxy core more slowly. 

In Figure flOl the orbit is along a polar path and moves in the +Z direction through the side of 
the caustic at (X, Y, Z, v x , v y , v z ) = (20.1 kpc, 0 kpc, 0 kpc, 0 km s _1 , 0 km s^ 1 , 200.0 km s^ 1 ), 
normal to the Galactic plane. The simulation time is 1.1 Gyr. The dwarf galaxy moves to almost 
the same position after 1.1 Gyr in both the caustic ring and log potentials. In this case the dwarf 
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Fig. 9.— Y vs. X (Galactocentric Cartesian coordinates) for N-body simulations of a dwarf galaxy 
that moves along an orbit in the +Y direction when it is at the position of the n = 2 caustic ring at 
(X, y)=(20.1, 0) kpc in a log (bottom left panel) and caustic (top left panel) ring halo potential. 
The simulations were run for 2 Gyr. The forward (dotted line) and backward (solid line) orbit is 
also plotted. The right panels show histograms of the number of bodies as a function of Galactic 
longitude l. 
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Fig. 10.— Z vs. X (Galactocentric Cartesian coordinates) for N-body simulations of a dwarf 
galaxy that moves along an orbit in the +Z direction when it is at the position of the n = 2 caustic 
ring at (X,Z) = ( 20.1, 0) kpc in a log (bottom left panel) and caustic ring (top left panel) halo 
potential. The simulations were run for 1.1 Gyr. The forward (dotted line) and backward (solid 
line) orbit is also plotted. The right panels show histograms of the number of bodies as a function 
of Galactic latitude b. 
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Fig. 11.— Y vs. X (Galactocentric Cartesian coordinates) for N-body simulations of a dwarf galaxy 
that moves along an orbit in the +Y direction when it is 1 kpc away from the n = 2 caustic ring 
at (X, Y)=(21.1, 0) kpc (top left panel) or (X,Y)=( 19.1, 0) kpc (bottom left panel) in a caustic 
ring halo potential. The simulations were run for 2 Gyr. The forward (dotted line) and backward 
(solid line) orbit is also plotted. The right panels show histograms of the number of bodies as a 
function of Galactocentric distance. Note that there is an excess of of debris at distances greater 
than 20 kpc in the top right histogram and an excess of of debris at distances less than 20 kpc in 
the bottom right histogram. 
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disrupts more slowly in log halo; the central peak in a caustic halo has about 20% of the total 
number of bodies and the log halo has about 40% of the total. The caustic ring may have the effect 
of stripping more bodies from the dwarf galaxy in its passage through the disk. 

In Figure fm the orbit is designed to move in the +Y direction of the Galactic plane just inside 
and outside of the caustic at ( X , Y, Z. v x , v y , v z ) = (19.1 kpc, 0 kpc, 0 kpc, 0 km s -1 , 200.0 
km s _1 , 0 km s^ 1 ) and ( X , Y, Z, v x , v y , v z ) = (21.1 kpc, 0 kpc, 0 kpc, 0 km s , 200.0 km s _1 , 
0 km s _1 ). These simulations were run for 2 Gyr. Here an interesting result is found. If the orbit 
is started just inside the n = 2 caustic ring, most of the stars are disrupted interior to the rings 
with distances less that 20 kpc. If we start just outside the ring the tidal debris is concentrated 
greater then 20 kpc just outside the ring. 


In summary, we find that the caustic rings do have subtle influences on the formation of tidal 
tails, but the differences are not dramatic. The study in this section is motivated by the proposal 
that the Monoceros Ring, a ring of stars about 20 kpc from the G alacti c center, is associate d 
with a dark matter caustic ring at 20 kpc (jNataraian &; Sikivid 120071 . citing iNewberg et all 120020 . 
Certainly we designed the orbits in Figures l9lfTTl to move along a caustic ring to study its effects, 
but rea listically dwarf g alaxies would not necessarily fall in to the Milky Way with such prescribed 
orbits. Xu et al. (120151 ) found that the observed stellar ring structures in the disk are associated 
with radial waves with peaks above and below the Galactic plane. Xu et al, (20151) find an upward 
peak or overdensity of stars at 10 kpc, a downward peak at 12-14 kpc, then an upward peak at 
16-18 kpc, and a downward peak at 20-24 kpc. These peaks are coincidentally near the positions 
of the predicted caustic rings at 10, 13, and 20 kpc. 


7. Summary and Conclusion 


In thi s paper we have presented the first study of the caustic ring model of the Milky Way 
halo ([Duff y & Sikiv ie 20081 ), or “Sikive halo”, in relation to astronomical observations. We show 
that the caustic ring model is consistent with observations of stars in the Milky Way galaxy. 


We presented the formalism for calculating the gravitational acceleration of a caustic ring halo. 
Our open source code for a caustic ring halo is available in the NEMO Stellar Dynamics Toolbox 
and the Milkyway@horne client repository. We use this code to show that the caustic ring halo 
reproduces a roughly logarithmic halo, with large perturbations near the rings. The caustic ring 
halo has very little mass within 2 kpc of the Galactic center, resulting in a gravitational acceleration 
that points away from the Galactic center and towards the inner-most rings for distances less than 
2 kpc. This outward acceleration may not be realized in practice, due to hierarchical merging and 
the presence of baryons. 

We show that the caustic ring halo can reasonably match the known Galactic rotation curve. 
We were not able to confirm or rule out an association between the positions of the caustic rings and 
oscillations in the observed rotation curve, due to insufficient rotation curve data. The parameters 
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for a Galactic model that includes a Hernquist bulge, Miyamoto-Nagai disk, and caustic halo are 
listed in Table [2j The values for the disk mass, scale length, and scale height; and bulge mass 
and scale length that we found are comparable to those found in other studies which fit Galaxy 
component models to observations (IDehnen fc Binnevlll998l : IWidrow fe Dubinskill2005l : iMcMillan 


20 111 : iKafle et al .11 20 id ). 


We explore the effects of dark matter caustic rings on dwarf galaxy tidal disruption with N- 
body simulations. Simulations of the Sagittarius dwarf galaxy in a caustic ring halo potential with 
disk and bulge parameters that are tuned to match the Galactic rotation curve, match observations 
of the Sgr trailing tidal tails as far as 90 kpc from the Galactic center. Like the NFW halo, 
they are however unable to match the leading tidal tail. None of the caustic, NFW, or triaxial 
logarithmic halos are able to simultaneously match observations of the leading and trailing arms of 
the Sagittarius stream. 


We further show that simulations of dwarf galaxies that move through caustic rings are qual¬ 
itatively similar to those moving in a logarithmic halo. The dwarf galaxy disrupts more or less 
slowly depending on whether it moves on a planar or polar orbit in the caustic and logarithmic 
halo potentials. In simulations of dwarf galaxies that move just inside or outside of a caustic ring, 
the tidal debris is shown to be disrupted mostly to the interior or exterior of the ring respectively. 
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